Hyperstaticity and loops in frictional granular packings 
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Abstract. The hyperstatic nature of granular packings of perfectly rigid disks is analyzed algebraically and through numerical 
simulation. The elementary loops of grains emerge as a fundamental element in addressing hyperstaticity. Loops consisting 
of an odd number of grains behave differently than those with an even number. For odd loops, the latent stresses are exterior 
and are characterized by the sum of frictional forces around each loop. For even loops, the latent stresses are interior and are 
characterized by the alternating sum of frictional forces around each loop. The statistics of these two types of loop sums are 
found to be Gibbsian with a "temperature" that is linear with the friction coefficient fl when fl < I. 
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ANALYSIS 

Every 2D granular packing is topologically isomorphic 
to a polyhedron. The mapping replaces each grain by 
a vertex, each contact between grains by an edge, and 
each of the granular packing's pore spaces (the regions 
surrounded by an elementary loop of grains) by a face of 
the polyhedron. A rigid boundary surrounding a granular 
packing is simply another node of the polyhedron located 
on the far side of a sphere. Thus, the polyhedron is the 
polygonization of a sphere and is convex, and Euler's 
formula implies 

k = 8 + ?.-e (1) 

where k is the number of contacts, g is the number of 
grains, X is the number of loops and e = 1 . The value of 
e is one less than the Euler Characteristic x for a convex 
polyhedron because the boundary has not been counted 
as a grain, although it is counted as a vertex in Euler's 
formula. With periodic boundaries around the packing, 
the resulting polyhedron is the polygonization of a torus 
and Euler's formula still applies with e — X —0. 

The average number of contacts per grain (the average 
coordination number) (Z), and the average number of 
grain per loop (G), are 



as, 
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Some manipulation of these along with Eq.[T]obtains, 
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Dropping the last term introduces negligible error for 
packings with rigid boundaries if ^ >> e. 

The stability equations for a quasi-static granular 
packings of round, frictional, 2D disks may be written 



£/if = 



(4) 



where /„ is a A'-dimensional vector consisting of all the 
contact normal forces in the packing, f the contact tan- 
gential forces, and w a 2§-dimensional vector consisting 
of the X and y components of the body forces of all the 
grains in the packing, typically having elements in the 
top half (representing the x components) and mg in the 
bottom half (representing the y components), where m 
is the grain mass and g is gravity. The top equation is 
for translational stability and the bottom equation is for 
rotational stability. () is a ^-dimensional null-vector be- 
cause body forces do not induce torques on these disks. 
Because /„ does not affect rotation, we may treat the sec- 
ond equation separately as if the grains are translationally 
frozen as in the case of a random gear network, has 
dimensions g x k with g < k per Eq. [3] so the system of 
rotations is underspecified (hyperstatic). The system may 
be made isostatic by adding A — £ equations, one per loop 
of grains to within the value of e. (e is related to rotation 
of the boundary and packing as a whole.) 

For an elementary loop of G grains (see Fig. [T] where 
G ~ 4-) embedded in a granular packing, the rotational 
equations of the grains in that loop may be written as. 



s^f = F 



(5) 



where the elements of the G-dimensional f are the tan- 
gential forces on the interior contacts of the loop, and the 
elements of F are the exterior tangential forces forces on 
each grain. The determinant of this matrix. 



detj:^ = 



2 if G odd 
if G even 



(6) 




FIGURE 1. Example of an elementary loop. /, are interior 
forces. Fa are exterior forces. 



indicates that odd and even loops are fundamentally dif- 
ferent. If G is odd, then the system of forces can be 
solved immediately. If G is even, then the matrix is 
singular with a null space dimensionality of one, so one 
of its rows must be replaced by a vector that spans the 
null space. We find that vector by replacing any row with 
(ai,fl2, . . . jAg) and insisting that the determinant be non- 
zero, and this obtains, 

f(-l)V,7^0 (7) 

The most symmetric treatment is to set a„ = (—1)". This 
row of the new matrix calculates the alternating sum of 
the tangential forces around the loop, 

i3=f (-!)"/„ (8) 

n 

which we call the alternating loop sum. The value of 
j3 is specified in the corresponding row of F. A linear 
elasticity model of this loop can verify that j3 equals the 
locked-in stresses in the loop, which are independent of 
the external forces F„. 

We cannot define an alternating sum around a loop 
when G is odd, but we can calculate the non-alternating 
sum of interior tangential forces around the loop, 

G 

CC=Y.fn (9) 

n=l 

for G odd or even. It turns out that 

f -jLtiFn if G Odd 
a={ (10) 
I jLtii~^)"Fn if G even 

These differences between odd and even loops reflect 
that even loops of gears are free to rotate and yet maintain 



a locked-in, interior stress, whereas odd loops of gears 
are frustrated and cannot turn and yet cannot maintain 
an interior stress. Analysis of multiple loops in a granu- 
lar packing shows that for each additional loop added to 
the packing, one more loop sum equation must be added 
to to make it square and non-singular: an alternating 
loop sum j3 if G is even, or a non-alternating loop sum a 
if G is odd. Thus, the tangential forces summed around 
loops (not the vector forces) are expected to be the fun- 
damental entity in understanding hyperstaticity. 

NUMERICAL SIMULATIONS 

To study the statistics of these loop sums, discrete (or dis- 
tinct) element modeling (DEM) has been performed. The 
DEM model comprises of a polydisperse assembly of 
frictional circular particles. Full details of this model are 
provided in 1 1 ] . This model has been employed to exam- 
ine the constitutive response of granular assemblies, in 
two dimensions, under a^ariety of compression and pen- 
etration tests (e.g. 1 1, 2, 3, 4]). The contact laws adopted 
are similar to other DEM simulations (e.g. [5]) which 
employ spring, dash-pot and friction slider to model in- 
teraction at contacts, as first proposed in [6]. The key dif- 
ference between this model and the classical DEM I6i] 
lies in the contact moment: the model employed here in- 
corporates a moment transfer, in accordance with the so- 
called modified distinct element method (MDEM) |6, 7]. 
The analysis of the prior section could be extended to 
include moment transfer with some loss of clarity. 

We performed a series of simulations involving a fric- 
tional granular packing in 2D. In each test, a granular 
assembly is created from 1681 circular particles, whose 
radii are chosen randomly from a uniform distribution 
between 0.1 mm and 0.15 mm. The particles are dropped 
into a box with dimensions of 10 mm x 10 mm, under 
gravity, with the coefficient of friction between particles 
initially set to /i = 10^. The assembly is then allowed to 
settle to a state where the kinetic energy is negligible. 

All walls are assumed to have the same material 
properties as the particles. Damping coefficients are as- 
signed according to the formulas: b" ~ 0.1 V'Wmm^", 
b' = 0.l^m,ni„k', ¥ = O.lRmin^/mmink'' where m,„,„ is 
the mass of the smaller particle. The discrete time step 
used in the numerical integration of the equations of mo- 
tion is assigned a value according to: Af = 0. 1 \Jm„,in/k". 

The theory outlined in the previous sections addresses 
packings with perfect rigidity, but numerical simulations 
always have some compression at the contacts, which is 
known to affect (Z) . To test the limits of the simulation 
and identify a range of parameters where the theory can 
be tested accurately, we performed fifteen simulations at 
a range of stiffnesses k'', k' and k" and Coulomb friction 
coefficient /i. A summary of all the stiffness constants 



TABLE 1. Spring stiffness values used for tfie normal 
and tangential contact force and the contact moment, ]<" , 
respectively. 



Test 



Normal k" 

(N/m) 



Tangential k' 

(N/m) 



Rotational f 

(N/m) 



1 


1.05 X 10^ 


3.50 X 


104 


3.50 X 


102 


2 


2.10 X 10^ 


7.01 X 


104 


7.01 X 


102 


3 


4.20 X 10^ 


1.40 X 


105 


1.40 X 


103 


4 


8.41 X 10^ 


2.80 X 


105 


2.80 X 


103 


5 


1.68 X 10^ 


5.60 X 


105 


5.60 X 


103 


6 


3.36 X 10^ 


1.12X 


10^ 


1.12X 


10-* 


7 


6.72 X 10^ 


2.24 X 


10^ 


2.24 X 


10^ 


8 


1.34 X 10^ 


4.48 X 


10^ 


4.48 X 


104 


9 


2.69 X 10^ 


8.97 X 


lO*" 


8.97 X 


10^ 


10 


5.38 X 10^ 


1.79 X 


10^ 


1.79 X 


105 


11 


1.08 X 10^ 


3.59 X 


10^ 


3.59 X 


105 


12 


2.15 X 10^ 


7.17 X 


10^ 


7.17X 


105 


13 


8.60 X 10** 


2.87 X 


10** 


2.87 X 


10^ 


14 


3.44 X lO** 


1.15 X 


lO"* 


1.15 X 


10^ 


15 


1.38 X 10'° 


4.59 X 


10^ 


4.59 X 


10^ 




FIGURE 2. (Colour online) (Z) as a function of Coulomb 
coefficient /i for several contact stiffnesses. Additional values 
of ¥ , 1<! and fc" were tested at /i = 10, demonstrating that 
(Z) 3.0, the isostatic value for perfect rigidity and infinite 
friction coefficient. 



used in each test is presented in TablefT] 

The first three tests (1-3) are performed as follows. 
In each test, we lowered the value for /i and let the 
particles settle again to a negligible kinetic energy. Again 
the value for /i is lowered and the system is left to settle. 
This process is repeated until pi is reduced to a value of 
10~^. The same process is used for each test, using the 
same values of pL, as selected in test 1. Note the amount 
by which jx is decreased from 10'' to 10^^ is not uniform. 
Previous simulations showed that significant changes did 
not occur until /i is less than one. The results for these 
tests are shown in Fig.|2] 

For very small values of the coefficient of friction, 
10^^ to 10^^, we observe a near constant value for 
the average coordination number. This value decreases 
with increasing particle rigidity. The (Z) then decreases 
rapidly around p, = 0.025 to ^ = 1, before saturating 
again to a near constant value. To ensure the trends are 
reproducible, we repeated the test for a k" value that is 
two orders of magnitude higher than that used in Test 1 . 
As simulation times proved prohibitively long for very 
large values of k", test 7 was run only from /i = 10^ down 
to ;U = 10-"^. 

As shown in Fig. IH we also performed an additional 
eleven tests (4-6, 8-15) for very high values of fc" to 
determine the limiting value for (Z) for perfect rigidity 
and infinite friction coefficient. Fig.[3]shows a plot of (Z) 
versus fc" for li — 10. As the normal stiffness coefficient 
is increased, (Z) approaches the isostatic limit of 3. The 
plunge in (Z) for A;" > 10^ can be attributed to an increase 
in the number of rattlers. 

As shown in Fig. [H we find that the behavior of 
(Z)max — (Z) is a power-law with jX. 
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FIGURE 3. (Z) as a function of stiffness yf for /x = 10 data 
points shown in Fig.|2] (Z) approaches the isostatic limit of 3.0 
as k" oo, however, the sudden plunge for k" > is probably 
due to a sudden growth in the number of rattlers. 
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FIGURE 4. (Z)n]ax — (Z) behaves as a power law of fl with 
exponent 0.30, and transitions to a plateau near fl = 0.1. The 
dashed lines are a guide to the eye. 
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FIGURE 5. Average value of loop sums (o!)o(](j (solid dots), 
(of)even (X's), and average value of alternating loop sums 
(jS)even (open circles) versus friction coefficient jl. Solid line 
is a power law in /x with unity exponent. Dashed line is a guide 
to the eye. 

ELEMENTARY LOOPS 

The problem of finding elementary loops is a well- 
studied problem in graph theory. Note that what we call 
loops in granular packings are called cycles in graph the- 
ory, not loops which are something different. The prob- 
lem faced is to find the minimum cycle basis of a graph. 
This is a set of cycles such that there is a minimal number 
of edges in each cycle and these cycles "combined" form 
a basis of the cycle space of the graph. We have chosen to 
implement an algorithm discussed in JK]. The algorithm 
consists of four steps: (1) find the shortest paths between 
every pair of vertices; (2) generate cycles using the paths 
found; (3) sort all cycles by length; (4) find all linearly 
independent cycles. 

Figure |5] shows the values of a averaged over all even 
loops (o;)even, and averaged over all odd loops (oc)odd, 
and the values of j3 averaged over all even loops (/3)even 
for all the simulations at one value of k" . 

In all cases the distribution of loop sums appears to 
be an exponential decay with the decay constant equal to 
the inverse of the average value of the loop sums. In other 
words, it appears to be a Gibbs distribution. An example 
for /I ~ 10^^ is shown in Fig.|6l To test how closely these 
distributions follow a pure exponential decay, the least- 
squares difference is calculated and summed over all 
bins in the distribution. In Fig.|7]7?^ is plotted for each 
of the three types of loop sums and for packings at each 
value of /i. It is found that the distributions become pure 
Gibbsian for /X < 1, which is the same place where (a) 
and (j3) become proportional to (cf. Fig. |5]l. /i (or 
a linear function of it) may be interpreted as the latent 
stress temperature in a packing below that limit. 
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FIGURE 6. Normalized distribution of loop sums over all 
odd loops in a packing with il = 10^^. Dashed line is an 
exponential decay. 
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FIGURE 7. F?' (sum of squares) norm to quantify how 
closely the distribution of (solid dots), aeven (X's), and 
jSeven (open circles) follows the Gibbs distribution for different 
values of /x. Dashed line is the expectation value for a perfect 
Gibbs distribution sampled as herein. 
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